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A hs tract- A cquisition  and  ultimate  tracking  of  a  weak  GPS 
signal  faces  several  technical  challenges,  notably,  possible  data 
bit  sign  reversal  every  20  ms  and  tolerable  frequency  error 
inversely  proportional  to  the  integration  interval.  Brute  force 
search  over  all  possible  combinations  of  the  unknown  values  is 
computationally  prohibitive.  Assisted  GPS  relying  on  external 
infrastructure  for  timing,  data  bit,  and  frequency  error 
information  is  costly.  Coherent  techniques  such  as  the  block 
accumulating  coherent  integration  over  extended  interval 
(BACIX)  have  recently  been  proposed  to  increase  coherent 
integration.  Although  efficient,  such  coherent  methods  may  still 
be  too  expensive  except  for  high-end  receivers  and  may  not 
maintain  the  SNR  performance  when  there  are  large  frequency 
changes  over  the  extended  integration  interval. 

In  this  paper,  we  set  forth  a  novel  method  that  utilizes  the 
semi-coherent  scheme  for  post-correlation  integration.  It  is 
named  as  block  accumulating  semi-coherent  integration  of 
correlations  (BASIC)  and  can  be  viewed  as  an  extension  of  the 
BACIX  algorithm.  Although  less  sensitive  than  coherent 
integration,  semi-coherent  integration  based  on  inter-block 
conjugate  products  is  computationally  more  efficient.  In 
addition,  it  can  handle  large  frequency  changes. 

The  BASIC  algorithm  is  first  formulated  in  the  paper. 
Computer  simulation  results  are  then  presented  to  illustrate  the 
operation  and  performance  of  the  BASIC  algorithm  for  joint 
estimation  of  the  initial  frequency,  chirping  rate  (rate  of  change 
in  frequency),  bit  sync,  and  bit  sign  pattern. 

I.  INTRODUCTION 

A  Global  Positioning  System  (GPS)  receiver  needs  to 
measure  the  distances  to  at  least  four  GPS  satellites  at  known 
orbital  locations  in  order  to  obtain  an  estimate  of  its  own 
position  plus  its  clock  error.  These  distances,  also  called 
pseudoranges,  are  obtained  from  the  measured  propagation 
delay  times  of  the  radio  signals  broadcast  by  the  orbiting 
satellites  to  the  GPS  receiver,  multiplied  by  the  speed  of  light. 
The  propagation  delay  time  is  measured  as  the  difference 
between  the  time  a  particular  event  of  the  radio  signal  arrives 
at  the  receiver  antenna  and  the  time  that  same  event  leaves 
the  satellite  antenna.  The  time  of  arrival  is  tagged  with  the 
receiver  clock,  which  may  be  biased  and  be  subject  to  large 
clock  drift,  particularly  in  inexpensive  commercial  receivers. 
The  time  of  transmission,  as  read  from  ultra-precise  atomic 
clocks  onboard  GPS  satellites,  is  encoded  in  the  navigation 
message.  Also  modulated  on  the  same  radio  signals  are  the 
ephemeris  data  which  must  be  demodulated  from  the  data  bits 
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to  reconstruct  the  satellite  orbital  positions  at  the  time  of 
transmission  [2,  3,  4,  5]. 

The  received  GPS  signals  must  attain  a  minimum  power 
level  to  ensure  a  GPS  solution,  which  is  routinely  met  when 
there  is  a  clear  line-of-sight  (LOS)  view  from  the  receiver 
antenna  to  GPS  satellites  above  the  horizon.  However,  when 
the  line-of-sight  view  between  the  receiver  and  a  satellite  is 
obstructed  (e.g.,  due  to  foliage,  mountains,  buildings,  or  other 
structures)  the  GPS  signal  strength  may  be  severely 
attenuated,  leading  to  position  fix  with  poor  accuracy  and 
even  breakdown  of  tracking  loops  inside  a  regular  receiver.  It 
is  therefore  desired  to  improve  GPS  receiver  sensitivity  to 
operate  on  GPS  signals  of  very  low  power  level  (or  weak 
signals)  to  satisfy  the  requirements  of  location-based  mobile 
e-commerce  and  emergency  call  location  (E911). 

One  approach  to  acquiring  weak  GPS  signals  is  assisted 
GPS  (AGPS).  The  AGPS  approach  relies  upon  a  wireless 
data  link  (or  other  means)  to  distribute,  in  real  time,  such 
information  as  time,  frequency,  navigation  data  bits,  satellite 
ephemeredes,  and  approximate  position  as  well  as  differential 
corrections  to  special  GPS  receivers  equipped  with  a  network 
modem  so  as  to  reduce  the  uncertainty  search  space,  to  help 
lock  onto  signals,  and  to  assist  navigation  solution.  This 
approach,  however,  comes  with  a  heavy  price  associated  with 
installing  and  maintaining  the  wireless  aiding  infrastructure 
and  services  required  to  provide  the  coverage. 

To  enable  weak  GPS  signal  acquisition,  one  known 
technique  at  the  receiver  end  is  to  extend  the  signal 
integration  time.  Coherent  integration  is  more  gainful  than 
the  non-coherent  counterpart  [11].  However,  when  the 
coherent  integration  interval  extends  beyond,  say,  20  ms, 
which  is  the  duration  of  one  data  bit  for  GPS  C/A-code,  the 
sign  reversal  of  navigation  data  bits  becomes  destructive  if  it 
occurs  in  the  middle  of  a  long  coherent  integration.  As  a 
result,  an  additional  signal  parameter,  namely,  the  data  bit 
transition,  has  to  be  searched  in  addition  to  the  two  usual 
parameters,  code  phase  (time)  and  carrier  frequency,  beside 
the  ID  number  of  GPS  satellites.  This  in  fact  constitutes  a 
four  dimensional  search. 

With  nominally  twenty-four  active  GPS  satellites  in  orbit, 
the  maximum  number  of  visible  satellites  is  about  ten  for  a 
near-Earth  receiver.  For  GPS  C/A-code,  the  amount  of  search 
in  time  (code  phase)  is  actually  fixed.  With  a  search  step  of  a 
half  code  chip,  there  are  2046  steps  needed  to  cover  the  entire 
code  sequence  of  1023  chips.  Sequential  stepping  through  all 
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the  2046  code  phases,  use  of  2046  correlators  in  parallel,  or  a 
combination  of  sequential  and  parallel  techniques  are  among 
a  host  of  approaches  of  different  complexity  for  code  search 
and  acquisition  [2,  3,  4,  5], 

However,  the  amount  of  search  in  frequency  and  in  data  bit 
transition  increases  as  the  coherent  integration  prolongs.  The 
frequency  uncertainty  for  a  stationary  user  (due  to  the  relative 
motion  between  the  user  and  a  GPS  satellite  and  the  receiver 
local  clock  drift)  is  typically  ±5  kHz.  When  the  integration 
time  is  increased  from  the  C/A-code  epoch  of  1  ms  to,  say,  20 
ms,  the  amount  of  frequency  search  needs  to  be  increased  by 
20-fold.  Similarly,  the  number  of  possible  locations  of  data 
bit  transition  also  increases  as  the  coherent  integration 
interval  increases.  The  number  of  possible  combinations 
increases  exponentially.  There  are  also  second-order  effects 
associated  with  the  extended  coherent  integration.  This 
mainly  involves  changes  in  the  carrier  frequency  and  code 
chipping  rate  due  to  acceleration  and  instability  of  local 
clock,  which  may  become  significant  over  long  integration 
intervals. 

A  typical  search  process  for  a  GPS  receiver  employing  1 
ms  correlation  for  the  C/A-code  using  either  sequential  or 
parallel  correlators  involves  search  in  four  dimensions.  The 
first  search  is  to  look  for  a  visible  GPS  satellite  among  those 
active  satellites  in  orbit.  For  a  given  satellite,  parametric 
search  is  conducted  in  three  nested  loops.  The  outer 
parametric  loop  is  typically  the  frequency  search,  which  steps 
through  the  frequency  uncertainty  interval  500  Hz  per  step. 
The  selected  frequency  offset  (a  frequency  bin)  is  added  to 
the  nominal  frequency  to  control  the  carrier  numerically 
controlled  oscillator  (NCO)  whose  output  drives  the  in-phase 
and  quadrature  components  of  the  reference  carrier  for  down- 
conversation  to  baseband. 

The  middle  parametric  loop  is  the  code  phase  search  over 
the  entire  code  epoch  of  1023  chips,  a  half  chip  per  step  (a 
code  lag)  typically.  In  a  correlator-based  receiver,  the  code 
phase  search  may  be  done  sequentially  one  code  phase  at  a 
time  or  in  parallel  with  all  code  phases  searched  at  the  same 
time. 

The  inner  loop  accumulates  the  1  ms  correlations  for  each 
search  grid  point  in  terms  of  a  frequency  bin  and  a  code  lag 
over  a  certain  period  of  time  (a  dwell  time),  during  which  the 
signal  is  added  up  for  detection.  Under  normal  signal 
strength,  the  dwell  time  is  usually  chosen  to  be  several  ms 
wherein  the  1  ms  correlations  are  power-combined  (i.e.,  the 
non-coherent  integration)  with  the  navigation  data  bit 
transitions  squared  out,  thus  not  being  an  issue.  However,  for 
weak  GPS  signals,  the  dwell  time  may  extend  beyond  20  ms 
for  coherent  integration  where  the  navigation  data  bit 
transition  plays  a  critical  role. 

When  the  despreading  correlation  is  carried  out  every  1 
millisecond,  the  complex  correlation  is  available  at  the  rate  of 
1  kHz.  To  further  integrate,  there  will  be  twenty  1 
millisecond  correlations  (20  data  points)  per  a  data  bit.  When 
coherent  integration  is  beyond  two  data  bits,  the  number  of 
possible  bit  transitions  is  large  and  bit  transition  patterns 
become  complicated.  There  are,  consequently,  outstanding 


problems  with  known  techniques  for  weak  GPS  signal 
acquisition  and  tracking  in  a  standalone  high-sensitivity 
receivers. 

Recently,  a  new  approach  for  weak  GPS  signal  acquisition 
without  network  assistance  was  presented  in  [10],  which 
belongs  to  the  category  of  approaches  aimed  at  enhancing  the 
sensitivity  of  standalone  GPS  receivers  (or  unaided  GPS 
receivers).  The  approach,  named  the  “block  accumulating 
coherent  integration  over  extended  interval  -  BACIX”  is 
based  on  three  ideas.  The  first  idea  of  successive  sign  reversal 
of  1  ms  correlations  per  block  (i.e.,  20  ms)  enables  data  bit 
transition  detection  and  data  bit  sign  correction  within  a  data 
bit  interval.  This  is  in  shaip  contrast  to  conventional  methods 
that  only  process  half  blocks  at  a  time  (undesired  SNR  loss) 
or  are  stretched  over  two  blocks  (excessive  latency).  In 
addition,  it  allows  for  efficient  use  of  Fast  Fourier  Transform 
(FFT),  which  only  needs  to  be  calculated  once  and 
manipulated  simply  for  all  tentative  bit  alignments.  As  the 
second  idea,  the  FFT  is  implemented  as  a  bank  of  bandpass 
filters  to  integrate  blocks  of  correlations  in  a  coherent 
manner,  as  opposed  to  conventional  phase  rotation  that  can 
only  handle  one  bandpass  filter  at  a  time  and  requires 
excessive  frequency  search.  Finally,  the  third  idea  optimizes 
the  overall  computation  by  pruning  unlikely  branches  of  the 
search  paths  in  partial  sums.  Such  a  block  accumulating 
coherent  integration  boosts  the  signal  power  while  reducing 
noise  despite  unknown  data  bit  transitions  and  other 
variations  that  may  occur  during  the  extended  coherent 
integration  intervals. 

Another  new  method  for  weak  GPS  signal  acquisition 
without  network  assistance  was  presented  in  [12],  which  is 
particularly  suitable  for  high  dynamic  applications.  In  the 
method,  the  post-correlation  signal  is  modeled  as  a  chiip  or 
linear  frequency  modulated  (LFM)  signal  wherein  the  effect 
of  high  dynamics  is  accounted  for  by  the  chirping  rate,  albeit 
unknown.  The  Wigner-Ville  distribution  (WVD)  is  used  to 
represent  the  signal  energy  in  the  time-frequency  plane  while 
the  Hough  transform  (HT)  integrates  the  energy  belonging  to 
a  chirp  signal,  which  distributes  along  a  line  in  the  time- 
frequency  plane  [1],  Similarly,  the  Radon  transform  can  be 
used  in  place  of  the  Hough  transform  to  integrate  the  chirp 
energy  along  a  line  in  the  time-frequency  plane  [7].  However, 
the  combination  of  the  WVD  and  HT,  the  Wigner-Hough 
transform  (WHT),  or  the  combination  of  the  WVD  and  RT, 
the  Wigner-Radon  transform  (WRT),  cannot  be  applied 
directly  to  the  despread  correlations  because  of  the  presence 
of  unknown  data  bits.  The  new  method  is  a  block- 
implemented  WHT/WRT  using  the  intra-block  conjugates 
products  [12].  They  were  shown  to  be  able  to  jointly  estimate 
the  initial  frequency,  chirping  rate  (frequency  rate),  and  bit 
sync  but  not  the  bit  sign  pattern. 

In  this  paper,  a  novel  approach  for  weak  GPS  signal 
acquisition  without  network  assistance  is  set  forth.  The 
approach  is  named  the  “block  accumulating  semi-coherent 
integration  of  correlations  -  BASIC,”  which  can  be  viewed  as 
an  extension  of  the  BACIX  algorithm  in  two  aspects.  First,  it 
can  handle  large  frequency  changes.  Second,  it  is 
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computationally  more  efficient.  The  novel  BASIC  algorithm 
is  based  on  the  semi-coherent  integration  scheme  as  analyzed 
in  [11]  and  makes  use  of  block  conjugate  products  [12]. 

After  introducing  a  complex  model  for  post-correlation 
GPS  signals  in  Section  II,  we  describe  the  BASIC  algorithm 
using  the  inter-block  conjugate  products  in  Sections  III.  In 
Section  IV,  simulation  results  are  presented  to  illustrate  the 
operation  and  performance  of  the  BASIC  algorithm  for  joint 
estimation  of  the  initial  frequency,  chirping  rate  (frequency 
rate),  bit  sync,  and  bit  sign  pattern.  A  summary  is  provided  in 
Section  V  with  concluding  remarks  and  future  work. 

II.  COMPLEX  MODEL  FOR  POST-CORRELATION  GPS  SIGNALS 

Consider  a  scheme  for  acquiring  GPS  signals  via  search  in 
time  and  frequency.  Assume  that  the  frequency  search  step  is 
500  Hz,  the  code  search  step  is  half  a  chip,  and  the  integration 
interval  of  the  despreading  correlation  is  1  ms.  Such 
integration  can  tolerate  ±250  Hz  errors  in  frequency  and  ±14 
chips.  For  a  search  grid,  the  worst  code  uncertainty  is  a 
quarter  of  a  chip  duration  {At  =  TJ4)  and  the  worst  Doppler 
uncertainty  is  Af=  250  Hz.  As  a  result,  the  maximum  loss  in 
a  1  millisecond  coherent  integration  (7)  =  1  millisecond)  is 

Loss  =  201ogioR(zk)  ±  201ogio  [sin( jrAfT/)/ 7iAfT/\ 

=  201oglo(0.75)  ±  201ogio(0.9)  =  -3.4  dB  (1) 

where  R{- )  stands  for  the  correlation  function  with  an  ideally 
triangular  shape. 

The  resulting  complex  correlation  is  unique  for  each  GPS 
satellite  and  can  be  written  as: 

exp{y[2^(/0«7;  +  an2Ts2)  +  ^0]}  ±  wn  (2) 

where  x„  is  the  complex  post-correlation  signal  with  the 
sample  index  n,  A„  is  the  amplitude,  (/>„  is  the  initial  phase,  Ts 
=  1  ms  is  the  integration  interval  with  the  sampling  rate  of  1 
kHz,  b„  =  ±1  is  the  unknown  data  bit,  fo  is  the  initial 
frequency,  a  is  the  chirping  rate  (sometimes  a/2),  and  w„  is 
an  uncorrelated  noise  with  zero  mean  and  unity  variance. 

For  GPS  C/A-codes,  the  baud  rate  is  50  Hz  and  each  data 
bit  thus  has  a  periodicity  of  20  ms.  There  will  be  20  samples 
per  data  bit  of  the  same  sign,  but  it  may  change  sign  from  one 
data  bit  to  the  next.  Knowing  which  sample  is  the  start  of  a 
data  bit  is  called  the  bit  sync.  Without  knowing  the  bit  sync,  a 
perfect  correction  may  be  destroyed  if  there  is  a  data  bit  sign 
reversal  in  the  middle  of  integration.  It  is  equally  important  to 
determine  the  sign  of  data  bits,  from  which  the  navigation 
messages  can  be  demodulated. 

The  instantaneous  frequency  in  (2)  represents  the 
frequency  error  between  the  incoming  signal  and  the  local 
carrier  replica,  which  can  be  written  as: 

Af„=f0nTs+an2Ts 2  (3) 

For  the  local  carrier  replica  fixed  at  a  search  frequency,  the 
complex  model  in  (2)  is  useful  as  long  as  the  frequency  error, 
albeit  being  time -varying,  remains  within  ±250  Hz. 
Otherwise,  the  search  switches  to  another  search  frequency. 


Over  a  time  interval  of  7],  the  maximum  chirp  rate  for  a  fixed 
carrier  replica  is: 


a  =  -signify) 


500- 1  /o  I 

T, 


(4) 


where  sign (•)  is  a  sign  function,  so  that  the  change  in 
frequency  is  less  than  500  Hz  in  one  second.  The  maximum 
chirp  rate  is  therefore  a  =  500  Hz/s.  At  LI  =  1575.42  MHz, 
the  frequency  change  over  one  second  is  related  to  the 
underlying  acceleration  (in  g  =  9.8  m/s2)  by  the  factor  of 
51.49  Hz/g/s. 

The  signal  amplitude  A„  includes  a  factor  7?(Ar)  related  to 
timing  error  A r  and  a  factor  sin(/r7(  A/')/ VTA/'  related  to 


frequency  error  A f  [2,  3,  4,  5],  As  a  result,  the  signal 
amplitude  A„  is  also  subject  to  variations  due  to  changes  in 
A f.  Since  the  change  is  rather  small,  the  amplitude  is  assumed 
to  be  constant  for  simplicity  in  this  paper. 

For  weak  signal  acquisition  and  other  applications,  it  is 
desired  to  further  integrate  the  1  ms  correlations,  xn,  over  a 
rather  long  period  of  time,  7).  Because  of  the  unknown  data 
bits  b„  and  varying  frequency  errors  A f  it  is  counter¬ 
productive  to  sum  complex  correlations  x„  directly.  Squaring 
can  move  both  the  data  bit  and  the  unknown  frequency  error 
and  allows  for  direct  summation.  However,  this  non-coherent 
integration  also  squares  noise  [11]  and  further  strips  off  all 
valuable  information  about  the  signal  that  is  needed  for 
subsequent  processing  in  a  GPS  receiver.  In  contrast,  a  semi- 
coherent  integration  scheme  is  used  to  implement  the  BASIC 
algorithm  based  on  inter-block  conjugate  products  as 
described  below. 


III.  BASIC  WITH  INTER-BLOCK  CONJUGATE  PRODUCTS 

Referring  to  Fig.  2,  consider  two  consecutive  blocks 
denoted  by: 

yt  =  [xki,  i  =  l,...,M]  =  [x{t)  t  =  {M{k -\)  +  i)Ts]  (5a) 

lk+l  =l>(r+i)/>  i  =  l,...,M]  =  [x(t  +  T),  t  =  {M{k-\)  +  i)Ts\  (5b) 

where  k  =  1,  ...,  K- 1  is  the  block  index,  v=  MTS  =  20  ms  is 
the  fixed  delay  between  two  blocks  with  for  M  =  20  and  Ts  = 
1  ms. 

Construct  a  (K-l)xM  matrix  of  inter-block  conjugate 
products  as: 


Z  =  [zki,  i  =  *  =  1,...,A:-1] 

II 

rcT 

'-t. 

II 

K> 

1 

+ 

(6a) 

zu  =  •••>  M,  k-  1,  ...,  K-\ 

(6b) 

z(t )  =  x{t  +  v)x\t )  , 

t  =  ( M(k-l)Ts ,  i  =  1,  . . .,  M,  k  =  1,  . . .,  K-\ 

(6c) 

With  the  chirp  signal  modeled  in  (2),  an 
conjugate  product  of  (6)  can  be  further  written  as: 

inter-block 

z(t)  =  x(t  +  t)x’  (t)  =  b(t)b(t  ±  r)zl  V2*(/»rW+2as) 

(7a) 
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Fig.  1.  Time  Diagram  for  Inter-Block  Conjugate  Products 
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The  complex  exponential  in  (7)  does  not  contain  any  t2 
term,  which  is  eliminated  by  the  delayed  conjugate 
products.  Putting  aside  the  unknown  sign  of  data  bit 
products  bkbk+l,  the  application  of  FFT  on  z/„  over  k  will 
produce  a  peak  at  2  ar  for  each  i,  from  which  the  chirp  rate 
can  be  estimated.  The  complex  peaks  can  be  FFT-analyzed 
again  over  i  to  further  strengthen  the  signal.  This  in  fact 
involves  a  2D  FFT  [6,  9]. 

Another  way  to  strengthen  the  signal  is  to  sum  up  all  the 
terms  pertaining  to  a  block  k  ( i.e .,  straight  summation  over  i 
for  each  k).  This  will  incur  a  small  loss  but  reduce 
computation  as  compared  to  the  use  of  FFT.  The  averaged 
term  can  be  written  as: 


z 


k 


=  btbk+1A- 


sm(2x;aM2T2) 

sm(27iaMT2) 


ej2tia(M+l)T;  gj2x(f„T+ar-  +2aT2M2  (*-!)) 


~bkbk+lMA2e 


j2tr(f„r+aT2  +2aT;M2  [k-l)) 


(8a) 

(8b) 


requires  determination  of  data  bit  signs  in  addition  to  the 
location  of  data  bit  transitions  (bit  sync).  One  simple 
approach  to  bit  sync  is  to  maintain  twenty  sums  of  blocks. 
The  start  of  blocks  of  a  sum  corresponds  to  a  possible  bit 
transition  location  within  20  ms  and  the  block  in  one  sum  is 
delayed  by  one  sample  (1  ms)  relative  to  the  preceding  sum 
as  shown  in  Fig.  2. 

The  klh  block  in  the  bth  sum  is  therefore  constructed  as: 

y_k  =  [*20(Jt-l)+i+6-l  >  f  =  lv>20]  =  i  = 

b  ^1,  ...,  20,  k=  1,  ...,K  (12) 

The  averaging  in  (8)  is  applied  to  each  and  every  block  of 
the  twenty  partial  sums.  These  averaged  blocks  are  then 
accumulated  with  bit  signs  rectified  using  the  method 
described  below.  At  this  point,  denote  the  twenty  sums  over 
K  blocks  by  Sb(vb,  n  =  1,  ...,  K).  The  partial  sum  that 

t-k 

produces  the  maximum  value  provides  an  estimate  of  the 
data  bit  transition  (i.e.,  the  tentative  bit  sync)  as: 

b  =  argmaxS,6(/,...,y\...V  )  (13) 

b  — 1  — n 


Then  applying  the  FFT  to  Zj  over  k,  without  considering 
the  data  bits,  produces: 

Zv=fk{[(zk,  k  =  l,2,...,K),  0^]}  (9) 

where  k  is  the  factor  of  K- 1  zeros  padded  to  z.  prior  to  the 
FFT. 

The  peak  location  of  Z,,  within  ±(k+\)K/2  is  obtained  as: 

/  =  arg  max  |  Zv  (10) 

v 

The  peak  location  /  provides  an  estimate  of  the  quantity 
2  ar,  from  which  the  chirping  rate  a  can  be  estimated 
directly  as: 

a(z)  =  — — — ^ — -  (11) 

2T2M-(K-\)(k  +  Y) 

Different  from  the  intra-block  conjugate  products  where 
the  product  of  two  samples  within  a  block  aligned  to  the 
same  data  bit  is  always  1  [12],  the  inter-block  conjugate 
products,  however,  may  have  different  signs.  This  therefore 


The  method  to  determine  data  bit  signs  and  accumulate 
the  partial  sums  with  bit  signs  rectified  is  now  presented. 
For  the  aligned  sum,  i.e.,  the  particular  sum  among  the 
twenty  delayed  sums  that  happens  to  have  the  start  and  end 
of  its  blocks  coincide  with  data  bits,  the  average  signal 
model  in  (8)  is  repeated  below  with  simplified  notations  as: 

zk=5kBeR2**+M  (14) 

where  8k  =  bkbk+\,  B  =  MA2,  y=  2 aT2M2,  =  2 7if0r+  ar2- 
2aTs:M2).  It  is  interesting  to  note  that  for  non-aligned 
delayed  sums,  8k  is  no  longer  binary-valued  because  of 
possible  sum  of  samples  across  sign  reversal. 

We  therefore  start  with  a  sequence  of  samples  ^  k  =  1, 

2,  ...,  K- 1,  available  at  50  Hz.  Each  sample  interval  is  the 
same  as  a  data  bit  of  20  ms.  To  determine  the  sign  for  each 
sample  8k  while  estimating  the  unknown  frequency  y,  we 
conduct  four  steps.  The  first  step  consists  of  pre  and  post 
zero-padding  each  sample  so  as  to  maintain  its  timing  order 
in  the  original  sequence: 


4 


**=[Q  L  z[  ^.^,*=1,2,. ..,^-1  (15) 

The  second  step  is  to  apply  the  FFT  to  the  K- 1  vector  hk 
expanded  from  a  single  sample  z  as: 


(16) 

where  /  is  the  index  of  the  frequency  bin  with  a  resolution 
inversely  proportional  to  the  data  length  ( K-\)MTS  (zero¬ 
padding  is  again  possible  if  a  refined  size  of  the  frequency 
bin  is  desired). 

The  third  step  is  to  obtain  the  twenty  partial  sums  up  to  k 
while  estimating  the  sign  Sk.  Assume  that  we  already  have 
estimated  the  bit  signs  up  to  k- 1  for  the  bth  partial  sum, 
denoted  by  #*(/),  and  obtained  the  corresponding  partial 

sums  up  to  k- 1 ,  denoted  by  (/) ,  for  /  =  1,2,  . . . ,  K- 1  and  b 

=  1,  2,  M.  Then,  the  partial  sums  and  bit  pattern  up  to  k 
are  obtained  by: 

sbk(i)=sU(!)  +  8bk(!)Hl(!)  (17a) 

Bbk(l)=[BUl),  %«)  j  (17b) 

#(0  =  argmax  |  #_,(/)  +  <£//*(/)  I’ 

■»,~U| 

/=  1,2,  1  and  b  =  1,2,  ...,M  (17c) 

In  the  fourth  step,  the  maximum  partial  sum  across  the 
frequency  bins  for  each  tentative  bit  transition  is  first 
obtained  as: 

lh  =  arg  max  \SbK_,(l)\^=  \,2,  ...,M  (18) 

Then,  the  maximum  partial  sum  across  all  tentative  bit 
transitions  is  obtained  as: 

/>  =  arg  max  |<S'£_1(/6)|  0^) 

which  provides  the  best  estimate  of  data  bit  transition,  thus 
achieving  the  bit  sync.  Finally,  the  best  frequency  bin  for 
the  maximum  partial  sum  at  the  estimated  bit  sync  is  given 
by: 

/  =  P  (20) 

and  the  corresponding  chirping  rate  a  can  be  estimated 
directly  from  (11). 

The  total  number  of  possible  sequences  to  maintain  is 
M(K- 1)  in  the  above  computation.  However,  it  is  possible  to 
trim  down  the  number  of  possible  bit  sequences  in 
frequency  bins  and/or  bit  transitions  after  few  samples  are 
integrated  by  discarding  the  smallest  partial  sums,  thus 
saving  computation  and  memory. 

Unique  to  the  inter-block  conjugate  products  is  the  fact 
that  the  sign  sequence  as  determined  above  is  a  differential 
bit  sequence  {8k  =  bkbk+\}.  The  original  bit  sequence  {bk} 
has  to  be  reconstructed  from  it.  To  illustrate,  consider  the 
example  listed  in  Table  1.  The  first  row  in  Table  1  shows 
the  original  bit  sequence  bk=  { 1,  -1,  -1,  1,-1,  1,  1,  ...}.  The 


second  row  is  the  differential  bit  sequence  is  bkbk+i  =  {-1,  1, 

-1,-1, -1,1,...}. 

The  inter-block  conjugate  products  allow  for  the 
estimation  of  the  differential  bit  sequence  8k  using  (17)  from 
the  average  samples  z  ■  However,  for  the  very  first  sample 

Z| ,  we  have  the  initial  conditions  s* (/)  =  ()•  Si  could  be 

either  1  or  -1.  If  the  first  selected  sign  matches  the  true 
differential  bit,  the  remaining  signs  will  reproduce  the 
underlying  differential  bit  sequence  exactly.  Otherwise,  it 
will  differ  as  being  the  complement. 

For  the  data  listed  in  Table  1,  the  first  sign  for  the 
estimated  differential  bit  is  fixed  at  1,  which  is  the  bolded 
number  in  the  third  row  of  Table  1 .  Since  it  differs  from  the 
first  true  sign,  the  estimated  differential  bits  (third  row)  are 
complementary  of  the  true  differential  bits  (second  row). 
Taking  complement  of  the  estimated  differential  bit 
sequence  (called  the  negative  logic)  in  the  fourth  row  turns 
out  to  be  the  same  as  the  true  differential  bits  (second  row). 

When  reconstructing  the  initial  bit  sequence  from  the 
estimated  differential  bit  sequence,  one  also  needs  to  select 
the  sign  for  the  very  first  bit.  Using  the  positive  and 
negative  logic  with  the  first  bit  fixed  as  “1”,  the 
reconstructed  bit  sequences  are  shown  in  the  fifth  and  sixth 
rows  of  Table  1,  respectively.  Clearly,  the  sequence  with 
negative  logic  is  correct  (the  sixth  row  vs.  first  row)  for  the 
example  shown  in  Table  1.  This  ambiguity  in  the  bit 
sequence  reconstruction  can  be  easily  solved  by  applying 
them  back  to  the  data  samples  and  the  correct  one  will  result 
in  the  same  sign  for  all  samples  and  the  resulting  sum  is 
bigger. 

For  the  signal  model  in  (14),  the  sign-rectified  sum  5*  (/) 
may  exhibit  two  peaks  of  equal  magnitude  at  frequency 
bins,  say,  /j  and  J_  for  two  differential  bit 
sequences  g^/j)  and  gj  q  ),  respectively.  When  applied  to 
the  samples  z  ,  the  correct  differential  bit  sequence  will 

render  the  signs  all  l’s  or  -l’s  whereas  the  incorrect 
sequence  makes  it  +1  and  -1  alternatively.  Assume  that  /j 
corresponds  to  the  correct  peak,  which  is  given  by 
S(l)exp{)27zjk)  =  exp{]27rixk}  because  g* _i(/j)  =  {3(1 ()  =  1 
or  -1  for  all  k}.  It  can  be  shown  that  the  second  peak 
location  is  related  to  the  correct  peak  location  by  /  =  /j  ±  14 
in  units  of  normalized  frequency.  Then,  the  peak  value  can 
be  evaluated  as: 

exp{j2xl2k}  =  exp{j2n(i{  ±  14 )k) 

=  exp {/2 k } exp { ±j 71k }  =  (+\)k+I exp{j2nlx  k}  (21) 

The  sign  term  of  (21)  is  cancelled  out  by  the  estimated 
differential  bit  g*  (/  )  =  ±1,  leading  to  the  same  value  as  the 
correct  peak.  Although  the  false  peak  parameters  are 
correctly  estimated,  they  can  be  distinguished  from  the  true 
signal  parameters  by  FFT-analyzing  the  de-chirped  de¬ 
modulated  signal  samples,  from  which  the  initial  frequency 
can  be  estimated  as: 
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Table  1.  Reconstruction  from  Differential  Bit  Sequence 


Initial  Bit  Sequence:  bk 

1-1-11-11  1 

Differential  Sequence:  b k  bk+\ 

/  -1  1  -1  -1  -1  1 

Estimated  Differential:  Sk 

1-11  1  1-1 

Complementary  (Negative)  Sk 

/  -1  1  -1  -1  -1  1 

Reconstructed  (Positive):  //’“ 

11  -1-1  -1  -1  1 

Reconstructed  (Negative):  i'”* 

1-1-11-11  1 

Table  2.  Bit  Sign  Determination 
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( f0,cc,b )  =  argmaxmax  |  FFT^b';  exp (-jl7din)xn,  n  =  1,..., iv]}  |  > 

li  >'  /o 

it  =1,2  and  i  =  1,2  (22) 

The  maximization  of  (22)  provides  the  global  peak  location 
indices  among  the  search  variables  i,  ih  and  the  FFT 
frequency  bin  f0.  The  frequency  bin  index  of  the  FFT  is  then 
translated  into  the  initial  frequency  estimate  f  ,  the  index  i  is 

translated  in  to  the  chirping  rate  a  via  i ,  and  the  logic  index 
ij  is  translated  into  the  estimated  bit  sequence  h . 

Clearly,  a  much  simpler  algorithm  results  if  the  signal  is  of 
low  dynamics,  i.e.,  a  «  0  Flz/s2.  The  application  of  FFT  over 
the  pre  and  post  zero-padded  samples  as  in  (16)  is  no  longer 
necessary.  This  will  considerably  save  computation  time.  The 
bit  sign  selection  remains  the  same  as  in  (17)  except  it  is  now 
only  applied  to  one  sequence  (i.e.,  I  =  0).  There  will  be  no 
ambiguity  due  to  double  peaks  but  one  still  needs  to 
distinguish  between  the  positive  and  negative  logic  used  in 
reconstructing  the  bit  sequence  from  the  estimated 
differential  bit  sequence.  This  can  be  easily  done  by  picking 
up  the  one  that  has  a  bigger  peak  in  the  FFT  applied  to  the 
corresponding  bit-wiped  samples. 

Between  the  two  extremes  (i.e.,  the  full  search  with  FFT  in 
(16)  and  /  =  0  for  a  «  0  as  described  above),  it  is  possible  to 
consider  a  small  number  of  chirping  rates  to  account  for 
medium  dynamics. 

In  summary,  the  post-correlation  integration  involves 
search  in  data  bit  transition,  data  bit  sign,  initial  frequency, 
and  chirping  rate.  Several  approximation  techniques  are 
presented  to  reduce  computation.  The  inter-block  conjugate 
products  transform  the  simultaneous  search  in  the  initial 


frequency  and  chirping  rate  into  two  one-dimensional 
searches  in  sequence.  Binary  bit  sign  search  is  done  with 
pinning  to  limit  the  tree  growth.  These  and  other  features 
make  the  BASIC  algorithm  computationally  viable  as  shown 
in  the  simulation  examples  below. 

IV.  SIMULATION  RESULTS  AND  ANALYSIS 

Simulation  examples  are  presented  below  to  illustrate  the 
operation  and  performance  of  the  method  using  inter-block 
conjugate  products  for  joint  estimation  of  signal  parameter 
estimation  and  bit  pattern  as  well  as  bit  sync.  In  the  first 
simulation  example,  we  also  consider  an  ideal  case  where 
there  is  no  noise  and  no  data  modulation.  The  integration 
interval  is  again  1  s  with  N=  1000  for  Ts  =  1  ms.  The  data  are 
divided  into  K  =  50  blocks,  each  with  M  =  20  samples.  The 
initial  phase  is  9=  32°,  the  initial  frequency  is^/o  =  -250  Hz, 
and  the  chirping  rate  is  all  =  250  Hz/s. 

Fig.  3  shows  the  2D  FFT  applied  to  the  inter-block 
conjugate  products  Z  as  defined  in  (20).  The  peak  value  is 
948.23  vs.  the  full  value  of  20x49  =  980  and  the  peak 
location  along  k  provides  an  estimate  of  the  chirping  rate, 
which  is  245  Hz/s  without  interpolation  vs.  the  true  value  of 
all  =  250  Hz/s.  The  apparent  errors  in  the  peak  value  and 
chirping  rate  are  due  to  quantization  error.  Either  zero¬ 
padding  prior  to  taking  FFT  or  curve-fitting  around  the  peak 
can  reduce  such  quantization  errors. 

Fig.  4  shows  the  FFT  applied  to  the  average  samples  zk 

with  k=  7  as  in  (9).  The  FFT  length  is  therefore  (k+\)xK  = 
8x49  =  392.  The  peak  value  is  963  vs.  the  full  value  of  980 
and  the  peak  location  provides  an  estimate  of  the  chirping 
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rate  according  to  (11),  which  is  248.725  Hz/s  vs.  the  true 
value  of  a  =  250  Hz/s.  As  indicated  by  (8),  the  sine  function 
is  clearly  visible  from  Fig.  4,  which  introduces  a  loss  to  the 
peak  value  and  sidelobes.  The  latter  may  be  attenuated  to 
some  extent  when  a  windowing  function  is  used. 

In  the  second  simulation  example,  all  the  conditions  are 
kept  same  as  in  the  first  simulation  except  data  bits  are  now 
modulated  on  to  the  samples.  In  the  example,  the  first  bit 
transition  occurs  between  the  16lh  and  17th  samples.  In  other 
words,  the  16th  sample  is  the  end  of  the  first  incomplete  data 
bit  while  the  17th  sample  is  the  start  of  a  second  data  bit. 

For  the  partial  sum  that  has  the  correct  alignment.  Fig.  5 
shows  the  bit  sign-rectified  FFT  of  the  block-averaged 
samples,  obtained  using  (17).  The  length  of  FFT  is  (k+\)xK 
=  392  for  k=  7.  As  analyzed  above,  there  are  two  peaks  of 
equal  strength  of  963  at  dt  =  248.73  Hz/s  and  d2  =  -994.90 
Hz/s,  respectively,  which  are  separated  by  a  half  of  the  x-axis 
scale  in  Fig.  5,  which  is  1250  Hz/s. 

For  each  peak,  there  is  an  associated  pattern  of  differential 
bits  determined  with  the  first  sign  fixed  at  1 .  As  explained  in 
Table  lor  each  differential  bit  pattern,  two  possible  bit 
sequences  are  reconstructed  using  a  positive  logic  and  a 
negative  logic,  respectively.  Table  2  summarizes  the  results. 
The  first  row  lists  the  first  20  signs  of  the  original  bit 
sequence.  The  second  row  shows  the  corresponding 
differential  bit  sequence.  The  next  three  rows  list  the 
estimated  differential  bit  sequence  associated  with  the  first 
peak  of  Fig.  5  with  the  first  sign  fixed  at  1  (3rd  row),  the 
reconstructed  bit  sequence  using  the  positive  logic  (4lh  row), 
and  the  reconstructed  bit  sequence  using  the  negative  logic 
(5th  row).  The  last  three  rows  list  the  estimated  differential  bit 
sequence  associated  with  the  second  peak  of  Fig.  5  with  the 
first  sign  fixed  at  1  (6th  row),  the  reconstructed  bit  sequence 
using  the  positive  logic  (7th  row),  and  the  reconstructed  bit 
sequence  using  the  negative  logic  (8th  row). 

To  test  which  one  of  these  four  possibilities  is  true,  the 
original  signal  samples  are  dechirped  with  the  estimated 
chirping  rate  «  and  d2  and  demodulated  with  the  bit 


sequences.  The  FFT  is  then  applied  to  the  dechirped  and 
demodulated  samples  with  the  results  shown  in  Fig.  6. 
Clearly,  the  negative  logic  bit  sequence  associated  with  the 
first  peak  (6th  row  in  Table  2)  is  the  correct  estimate  (except 
for  a  complementary  sign  difference).  The  peak  location  at 
frequency  bin  751  provides  an  estimate  of  the  initial 
frequency  as  (751-1000)  =  -249  Hz  vs.  the  true  value  of /(l  =  - 
250  Hz.  The  false  peaks  have  50  Hz  harmonics  reflecting  the 
underlying  data  bits.  In  the  third  simulation  example,  all  the 
conditions  are  kept  same  as  in  the  second  simulation  with 
complex  noise  added  such  that  the  resulting  SNR  is  -3  dB. 
Under  the  Nyquist  rate,  the  signal  passes  through  a  low  pass 
filter  with  an  equivalent  bandwidth  of  ±//./2.  For  Ts  =  0.001  s, 
SNR  =  -3  dB  results  from  A  =  0.7,  which  also  corresponds  to 
C/N0  =  27  dB-Hz. 

Again,  for  the  partial  sum  that  has  the  correct  alignment, 
Fig.  7  shows  the  bit  sign-rectified  FFT  of  the  block-averaged 
samples.  Fig.  8  shows  the  FFT  of  dechirped  and  demodulated 
samples.  Compared  to  Figs.  5  and  6,  Figs.  7  and  8  show 
signal  peaks  lowered  by  a  factor  of  two,  consistent  with  SNR 
=  -3  dB.  However,  the  peaks  remain  at  the  same  locations, 
thus  providing  the  correct  estimates  of  chirping  rate,  initial 
frequency,  bit  pattern,  and  bit  sync. 

Fig.  9  shows  the  twenty  delayed  sums  for  the  case  without 
noise  over  an  interval  of  3  s.  The  peak  is  around  the  correct 
bit  transition  location  and  the  valley  is  about  10  samples 
away  from  the  peak  (partial  sums  starting  around  the  middle 
of  a  data  bit).  However,  the  peak  is  not  as  sharp  as  the 
method  using  the  intra-block  conjugate  products  [12].  This 
may  be  due  to  the  fact  that  complex  exponentials  across  a 
block  (vs.  within  a  block)  prevent  samples  that  have  opposite 
signs  from  complete  cancellation.  The  flattened  peak  makes 
precise  bit  sync  difficult  over  a  short  data  span. 

Fig.  10  shows  the  twenty  delayed  sums  for  SNR  =  -3  dB 
over  a  1  s  interval.  The  cumulated  peak  corresponds  to  the 
true  bit  transition,  thus  achieving  bit  sync.  However,  we 
observed  cases  where  the  peak  location  is  with  two  samples 
(±2  ms)  of  the  true  bit  transition.  To  note,  in  practical 
systems,  it  may  take  thousands  of  bits  (in  seconds  to  minutes) 
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to  firmly  establish  bit  sync,  followed  by  continuous 
monitoring  of  bit  transitions  and  bit  re-sync  if  necessary. 

V.  CONCLUSIONS 

In  this  paper,  the  GPS  signal  after  the  1  ms  despreading 
integration  was  modeled  as  a  modulated  chirp  signal  where  in 
the  binary  modulation  represents  the  data  bits  of  navigation 
message  while  the  chirp  accounts  for  changes  in  residual 
Doppler  frequency. 

To  further  integrate  the  signal  over  a  rather  long  period  of 
time,  the  semi-coherent  integration  scheme  (BASIC)  was 
used  on  the  inter-block  conjugate  products  of  complex 
correlations.  As  the  simulation  results  showed,  it  allowed  for 
joint  estimation  of  the  signal  parameters  (i.e.,  the  initial 
frequency  and  chirping  rate)  and  data  bits  (bit  sync  and /  bit 
signs).  There  were  ambiguous  bit  patterns  and  it  could  be 
easily  resolved  with  an  additional  step  of  tests. 

The  block-integrating  methods  presented  in  this  paper  can 
handle  large  frequency  changes  caused  by  relative  motion 
and  clock  drift.  It  can  also  remove  data  bits  in  an  open-ended 
manner,  allowing  for  new  samples  to  be  added  on  in  the 
integration  process  as  they  become  available.  Indeed,  when 
the  BASIC  algorithm  of  the  present  paper  is  applied  in  the 
downstream  together  with  the  MAGIC  (i.e.,  the 
multiresolution  adjustable  grid  interpolated  correlator)  [8]  in 
the  upstream,  a  new  GPS  receiver  architecture  emerges  with 
open-loop  acquisition  and  estimation.  It  is  in  contrast  to  the 
conventional  receiver  architecture  that  uses  closed-open 
acquisition  and  tracking.  This  open-loop  acquisition  and 
estimation  architecture  is  particularly  suitable  for  software 
radio  receivers  with  snapshot  solutions.  Research  is  under 
way  to  further  advance  this  concept. 
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